function [matrix_ACC_train,matrix_PPV_train,matrix_NPV_train,matrix_TP_train,matrix_FP_train,matrix_TN_train,matrix_FN_train,activation_weeks_var1_train,activation_weeks_var2_train] = fun_train(vec_training_time,windows_training,vec_bg,var1_timeseries,vec_week_train_event,vec_alpha)


%% Start loop on alpha_tol and bg  values

% Loop on GT data for the selected county. Compute selected windows in
% which there is an increase in the google trends variable above the
% chosen threshold OR there are a minimum of alpha_tol weeks of increasesin
% the data


% Start loop on digital data
for ind_alpha = 1:length(vec_alpha)
    tol_alpha = vec_alpha(ind_alpha);
    
    for ind_bg = 1:length(vec_bg)
        
        bg     = vec_bg(ind_bg); % set threshold for digital data
        count1 = 0;          % set a counter of events
        
        % Start loop on training windows to compute and store activation weeks for the given bg
        %threshold.
        for j1=1:length(windows_training)
            
            % For the threshold approach
            aux0 = var1_timeseries(windows_training{j1}) >= bg*max(var1_timeseries(vec_training_time));                      % pick up gt time series on the chosen window. Check which weeks are below or above-equal threshold
            aux1 = strfind(aux0',[0 1]) ;                                                                                   % indicates how many crossings (sequnce of [0 1] in aux0) happened: aux1 gives the list where the zeros of [0 1] are located.
            
            % For the derivative approach
            aux_windows_training     = windows_training{j1};                                                                    % For the derivative approach: store window with an auxiliary variable
            aux_alpha_vec_var1       = var1_timeseries(aux_windows_training(2:end))./var1_timeseries(aux_windows_training(1:end-1)); % compute the ratios of the 1-week lag time window over the current window. Build vector of x(t+1)/x(t) for the time series x.
            
            
            % If there is a crossing but not sufficiently increase in data
            if ~isempty(aux1) && length(find(aux_alpha_vec_var1>1))<tol_alpha                                                   % criteria for event detection: signal should have crossed the threshold at least at 1 point in the sliding window
                count1                                 = count1 + 1;                                                          % count as an event
                aux_select_windows1                    = windows_training{j1};                                                % store selected window
                select_windows1{ind_alpha,ind_bg,count1}             = aux_select_windows1    ;                                             % save window
                activation_weeks_var1_train{ind_alpha,ind_bg,count1} = aux_select_windows1(aux1(1)+1);                                      % save week where first crossing happened. This is the activation week.
                
                % If there is no crossing but sufficiently many increases in data
            else if isempty(aux1) && length(find(aux_alpha_vec_var1>1))>=tol_alpha
                    count1                        = count1 + 1;                                                         % count as an event
                    aux_select_windows1           = windows_training{j1};                                               % store selected window
                    select_windows1{ind_alpha,ind_bg,count1}    = aux_select_windows1;                                                % save window
                    
                    aux_alpha_vec_var1                    = find(aux_alpha_vec_var1>1);                                            % check index where an increase (alpha>1) happened
                    aux_activation_alpha_var1             = aux_alpha_vec_var1(tol_alpha);                                    % get ''tol_alpha''-th week where the increase started. Picks t if x(t+1)/x(t+1)>1
                    activation_weeks_var1_train{ind_alpha,ind_bg,count1} = aux_select_windows1(aux_activation_alpha_var1 + 1);               % The activation week is the next week after ''tol_alpha''-th week where the increase started. Picks t+1 if x(t+1)/x(t)>1
                    
                    %  If there is a crossing AND sufficiently many increases in data
                else if ~isempty(aux1) && length(find(aux_alpha_vec_var1>1))>=tol_alpha
                        
                        count1                         = count1 + 1;                                                         % count as an event
                        aux_select_windows1            = windows_training{j1};                                               % store selection window
                        select_windows1{ind_alpha,ind_bg,count1}  = aux_select_windows1;                                                % save window
                        activation_weeks_var1_cross    = aux_select_windows1(aux1(1)+1) ;                                    % save week where first crossing happened. This is the activation week.
                        
                        aux_alpha_vec_var1            = find(aux_alpha_vec_var1>1);                                           % check index where an increase (alpha>1) happened
                        aux_activation_alpha_var1     = aux_alpha_vec_var1(tol_alpha);                                   % get ''tol_alpha''-th week where the increase started. Picks t if x(t+1)/x(t+1)>1
                        activation_weeks_var1_alpha   = aux_select_windows1(aux_activation_alpha_var1 + 1);               % The activation week is the next week after ''tol_alpha''-th week where the increase started. Picks t+1 if x(t+1)/x(t+1)>1
                        
                        activation_weeks_var1_train{ind_alpha,ind_bg,count1} = min(activation_weeks_var1_cross,activation_weeks_var1_alpha);    % When two events happened in the chosen window, pick the one that happened first.
                        
                    end
                end
            end
            clear aux0 aux1 aux_alpha_vec aux_windows_training alpha_vec_var1
            
        end
        
        clear bg count1
        
    end
    
    clear tol_alpha
end



clear ind_alpha ind_bg j1 activation_weeks_var1_cross activation_weeks_var1_alpha ...
    aux_activation_alpha_var1 aux_activation_alpha_var1 aux_select_windows1 aux_alpha_vec_var1


if ~exist('select_windows1','var')
    matrix_ACC_train = [];
    matrix_PPV_train = [];
    matrix_NPV_train = [];
    matrix_TP_train = [];
    matrix_FP_train = [];
    matrix_TN_train = [];
    matrix_FN_train = [];
    activation_weeks_var1_train = [];
    activation_weeks_var2_train = [];
   return
end

%% define activation_weeks_var2_train as the cell with ''window_gap''entries of the event week
window_gap = length(windows_training{1})-1;

for j0=1:length(vec_week_train_event)
    
    train_week_event = vec_week_train_event(j0);
    
    if train_week_event>1
        
        aux_count2 = sort(1:min(window_gap,train_week_event-1),'descend');
        
        for count2 = 1:min(window_gap,train_week_event-1)
            aux_select_windows2{j0,aux_count2(count2)} = max(train_week_event-count2,1):train_week_event + window_gap - count2;
            aux_activation_weeks_var2_train{j0,count2}    = train_week_event;
        end
    end
    
    clear aux_count2 train_week_event
    
end

select_windows2            = reshape(aux_select_windows2',[1,j0*count2]);
activation_weeks_var2_train = reshape(aux_activation_weeks_var2_train',[1,j0*count2]);



%% COMPUTE ACC, PPV AND NPV

% Start loop on alpha_tol values
for ind_alpha = 1:size(select_windows1,1)
    
    % Start loop on bg threshold values and store ACC, PPV and NPV
    for ind_bg = 1:size(select_windows1,2)
        
        % Store info for confusion matrix
        num_TP  = 0; num_FP  = 0; num_TN  = 0; num_FN  = 0;
        
        % Redefine selected windows for digital data (var1) for scanning of TP, TN, FP, FN
        for ii1=1:size(select_windows1,3)
            select_windows1_new{ii1} = select_windows1{ind_alpha,ind_bg,ii1};
        end
        
        % Start Loop on windows: compute TP, TN, FP and FN
        for j3=1:length(windows_training)
            
            if  ~any(cellfun(@(x) isequal(x, windows_training{j3}), select_windows1_new)) &&  ~any(cellfun(@(x) isequal(x, windows_training{j3}), select_windows2))            % if there is no var1 AND no var2 activation at the given window, count as true negative
                num_TN = num_TN+1;
            else if ~any(cellfun(@(x) isequal(x, windows_training{j3}), select_windows1_new)) &&  any(cellfun(@(x) isequal(x, windows_training{j3}), select_windows2))         % if there is no Var1 activation but var2 activates, counts as false negative
                    num_FN = num_FN  + 1;
                else if any(cellfun(@(x) isequal(x, windows_training{j3}), select_windows1_new)) &&  ~any(cellfun(@(x) isequal(x, windows_training{j3}), select_windows2))     % if there is no Var2 activation but var1 activates, counts as false positive
                        num_FP = num_FP  + 1;
                    else if any(cellfun(@(x) isequal(x, windows_training{j3}), select_windows1_new)) &&  any(cellfun(@(x) isequal(x, windows_training{j3}), select_windows2)) % if both Var1 and Var2 activates, we ensure that var1 activates BEFORE var2 to count as True positive. Otherwise, we count as false negative. We continue as follows (se lines below)
                            
                      
                            % Start Loop on selected windows for var1
                            for jj1 = 1:length(select_windows1_new)
                                if isequal(windows_training{j3},select_windows1_new{jj1}) % check if the window is exactly a selected window of index iind_bg
                                    dummy_ind1(jj1) = 1;
                                else
                                    dummy_ind1(jj1) = 0;
                                end
                            end
                            
                            act_week1 = activation_weeks_var1_train{ind_alpha,ind_bg,dummy_ind1==1}; % find the activation week.
                            
                            % Start Loop on selected windows for var2
                            for jj2 = 1:length(select_windows2)
                                if isequal(windows_training{j3},select_windows2{jj2}) % check which if the window is exactly a selected window of index ii2
                                    dummy_ind2(jj2) = 1;
                                else
                                    dummy_ind2(jj2) = 0;
                                end
                            end
                            
                            act_week2  = activation_weeks_var2_train{dummy_ind2==1};% find the activation week.
                            clear dummy_ind1 dummy_ind2
                            
                   
                            % Check for TP or FN.
                            if act_week1 < act_week2
                                num_TP = num_TP+1;
                            else
                                num_FN = num_FN+1;
                            end
                            
                        end
                    end
                end
            end
            
            clear dummy_ind1 dummy_ind2 aux_var1 aux_var2  act_week1  act_week2
        end
        
        % Compute ACC, PPV and NPV
        ACC = (num_TP + num_TN)/(num_TP + num_TN + num_FP + num_FN); % accuracy
        PPV = num_TP/(num_TP+num_FP);                                % Positive predicted values
        NPV = num_TN/(num_TN + num_FN);                              % Negative predicted values
        
        % Store information in matrices
        matrix_ACC_train(ind_alpha,ind_bg) = ACC;
        matrix_PPV_train(ind_alpha,ind_bg) = PPV;
        matrix_NPV_train(ind_alpha,ind_bg) = NPV;
        matrix_TP_train(ind_alpha,ind_bg)  = num_TP;
        matrix_FP_train(ind_alpha,ind_bg)  = num_FP;
        matrix_TN_train(ind_alpha,ind_bg)  = num_TN;
        matrix_FN_train(ind_alpha,ind_bg)  = num_FN;
        
        clear   num_TP num_TN num_FP num_FN ACC PPV NPV  select_windows1_new select_windows2_new iind_bg ii2 j3 jj1 jj2
        
    end
    
end

end

